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ABSTRACT 


An  approach,  based  on  recent  work  by  Stern 
[STER  1979]  ,  is  described  for  obtaining  the  approxi¬ 
mate  transient  behavior  of  both  the  M/M/1  and  M(t)/M/1 
queues,  where  the  notation  M(t)  indicates  an  exponen¬ 
tial  arrival  process  with  time-varying  parameter  A(t). 
The  basic  technique  employs  an  M/M/l/K  approximation 
to  the  M/M/1  queue  to  obtain  a  spectral  representation 
of  the  time-dependent  behavior  for  which  the  eigen¬ 
values  and  eigenvectors  are  real. 

Following  a  general  survey  of  transient  anal¬ 
ysis  which  has  already  been  accomplished,  Stern's 
M/M/l/K  approximation  technique  is  examined  to  deter¬ 
mine  how  best  to  select  a  value  for  K  which  will  yield 
both  accurate  and  computationally  efficient  results. 

It  is  then  shown  how  the  approximation  technique  can 
be  extended  to  analyze  the  M(t)/M/1  queue  where  we 
assume  that  the  M(t)  arrival  process  can  be  approxi¬ 
mated  by  a  discrete ly  time-varying  Poisson  process. 

An  approximate  expression  for  the  departure 
process  of  the  M/M/1  queue  is  also  proposed  which  im¬ 
plies  that  for  an  M(t)/M/1  queue  whose  arrival  process 
is  discretely  time-varying,  so  too  the  departure  proc¬ 
ess  can  be  approximated  as  discretely  time-varying 
(albeit  with  a  different  time-varying  parameter). 

In  all  cases,  the  accuracy  and  validity  of 
these  results  are  examined  by  comparison  with  exact 
analytic  results,  simulation  or  alternative  discrete¬ 
time  approaches  (e.g.,  the  embedded  Markov  chain  tech¬ 
nique  of  Moore  |MOOR  1972]  (MOOR  1975]). 


i  ii 


1. 


INTRODUCTION 


The  concepts  of  stationarity ,  ergodicity,  and  conver¬ 
gence  are  fundamental  to  the  analysis  of  queues  in  equilibrium. 
The  basic  result  of  such  analysis  is  that  it  is  possible  to 
infer,  from  a  finite  number  of  observations  of  the  arrival  and 
departure  processes  of  a  queue,  the  stochastic  properties  of 
the  queue  over  all  time.  Processes  that  occur  in  real  life, 
however,  seldom  satisfy  these  conditions  and  consequently  do 
not  lend  themselves  to  steady  state  analysis. 

Despite  this  fact,  the  study  of  the  transient  behavior 
of  queues  has  been  largely  ignored.  The  reasons  for  this 
include  the  facts  that: 

1)  Very  few  transient  solutions  currently 
exist . 

2)  Those  that  do  exist  are  extremely  compli¬ 
cated  to  derive  and  manipulate. 

3)  Solutions  are  generally  limited  to  single 
queues . 

Transient  analysis  has,  moreover,  come  to  be  viewed 
in  two  somewhat  different  ways.  One  view  is  that  it  deals 
primarily  with  the  analysis  of  stationary  queues  approaching 
equilibrium.  In  this  context,  transient  queueing  analysis  is 
concerned  with  deriving  the  behavior  of  queues  as  they  converge 
to  equilibrium  as  well  as  the  speed  with  which  the  convergence 
occurs  (see,  for  example,  (JAIS  1960]  [TAKA  1961]  (FELL  1968] 

| BHAT  1968]  [ BS  1969]).  A  fundamental  concept  underlying  this 
version  of  transient  analysis  is  the  relaxation  time  (also 
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called  the  dominant  time  constant  [GR1F  1978]  [WSB  1975])  which 
is  a  measure  of  the  time  it  takes  a  queue  to  "achieve"  steady 
state.  Paradoxically,  the  efforts  put  into  this  form  of  tran¬ 
sient  analysis  have  been  generally  directed  at  justifying  the 
use  of  steady  state  solutions  by  showing  that  the  relaxation 
time  is  negligible. 

Another  view  of  transient  analysis  is  that  it  is  con¬ 
cerned  with  studying  the  behavior  of  queues  whose  input  and 
service  processes  are  time  dependent;  more  precisely,  the  pa¬ 
rameters  of  the  input  and  service  processes  are  held  to  be 
functions  of  time.  Such  queues  are  obviously  nonstationary 
and  will  never  achieve  equilibrium.  The  techniques  for  ana¬ 
lyzing  the  behavior  of  nonstationary  queues  fall  into  two  major 
categories:  those  which  treat  the  parameters  of  the  arrival 

and  service  processes  as  continuous  functions  of  time,  and 
those  which  treat  (or  approximate)  the  parameters  as  discrete 
(step)  functions  of  time. 

A  number  of  the  techniques  treating  the  parameters  as 
continuous  functions  of  time  rely  upon  the  Kolmogorov  difference- 
differential  equations  for  birth-death  processes.  Using  these 
equations,  Reuter  and  Ledermann  |RL  1953]  examined  the  M(t)/M(t)/1 
queue,  where  the  argument  t  indicates  that  the  process  has  a 
time  varying  parameter.  They  represented  the  queue  by  a  Markov 
process  with  an  enumerable  set  of  states,  and  then  provided 
sufficient  conditions  for  the  existence  of  a  unique  solution 
to  the  transition  probability  P^j^t),  the  conditional  proba¬ 
bility  that  the  system  is  in  state  k  at  time  t,  given  that  it 
was  in  state  i  at  time  0.  These  conditions  include  any  one  of 
the  following: 
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The  and  ^  are,  respectively,  the  birth  and  death  transition 
rates  associated  with  the  underlying  Markov  process. 


Takacs  [TAKA  1955]  investigated  the  virtual  waiting 
time  process  of  a  single  server  queue.  For  such  a  queue  with 
nonhomogeneous  Poisson  arrivals  (i.e.,  with  parameter  A(t)) 
and  general  service  time  distribution,  H(x),  he  was  able  to  derive 
the  integro-dif ferential  equation  for  the  waiting  time  dis¬ 
tribution,  F(t,x),  at  time  t: 


3F(t,x) 
3 1 


3F( t ,x ) 
3x 


-  A(t)F(t,x)  +  A ( t ) 


H(x-y )dyF(t ,y) 


He  then  derived  an  expression  for  $(t,s),  the  Laplace-Stielt jes 
transform  of  the  waiting  time: 


<J)(  t ,  s  )  =  e 


st-  [  l-4>(s)]A(t) 


1  -  s 


i 


F(u,0)  e 


■su+[  1  -ip ( s )  ]A(u) 


du 
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where  »|>(s)  is  the  Laplace-Stielt jes  transform  of  the  service 

time  distribution  and  A(t)  =  A.(u)  du.  Takacs  showed  that 

if  the  emptiness  probabilities  for  the  system  can  be  determined 
uniquely,  then  the  waiting  time  distribution  follows  directly. 

He  also  noted,  however,  that  deriving  the  emptiness  probabili¬ 
ties  is  usually  quite  difficult. 

Clark  [CLAR  1956]  attempted  to  solve  the  forward  Kol¬ 
mogorov  equations  for  the  M(t)/M(t)/1  queue  using  a  generating 
function  approach.  No  exact  representation  of  the  queue  length 
distribution  was  derived,  but  the  problem  was  "reduced"  to 
obtaining  the  solution  of  an  integral  equation.  Expressions 
were  derived,  however,  for  the  mean  and  variance  of  the  queue 
length  which  depend  on  the  emptiness  probabilities  for  the 
system.  Explicit  results  for  the  queue  length  were  obtained 
only  for  the  special  cases  where  A.(t)/p(t)  =  c,  some  constant. 

Another  continuous  time  approach  (proposed  by  Luchak 
I  LUCH  1956])  involves  the  queue  M(t)/E^/1,  where  E^  here  denotes 
a  weighted-sum  erlang  distribution  whose  integer  parameter  is 
a  non-negative  random  variable  Y.  Luchak  found  that  if  the 
arrival  process  parameter  can  be  expressed  as  a  polynomial  in 
time,  then  the  queue  length  distribution  can  be  formulated 
recursively;  only  in  selected  cases,  however,  can  a  closed 
form  solution  be  found.  Luchak  also  noted  that  periodic  ar¬ 
rival  rates  can  be  dealt  with  more  easily  than  most  since  the 
queue  length  distribution  need  only  be  found  for  the  initial 
period.  The  state  at  the  end  of  each  period  can  then  be  used 
as  the  initial  condition  for  the  next  period.  He  further  sug¬ 
gested  that  if  the  initial  state  of  a  queue  with  periodic  ar¬ 
rivals  was  taken  to  be  the  steady  state  solution  of  a  queue 
with  constant  arrival  rate  equal  to  the  average  rate  over  a 
period,  then  a  ’quasi-stationary"  steady  state  could  be  reached 
in  only  a  few  periods. 
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Keilson  and  Kooharian  (KK  1960)  examined  a  multi-di¬ 
mensional  phase  space  model  for  a  queue  with  time-dependent 
Poisson  arrivals  and  general  service  time  distribution,  and 
later  extended  their  analysis  to  general  time-dependent  arri¬ 
vals  [KK  1962). 

Hasofer  [HASO  1964]  [  HASO  1965],  using  the  results 
of  Reich  [REIC  1958]  [REIC  1959]  in  which  it  was  shown  that 
the  time  dependent  emptiness  probabilities  are  the  unique  so¬ 
lution  of  a  Volterra  equation  of  the  first  kind,  derived  expli¬ 
cit  results  for  the  M(t)/G/1  queue.  Assuming  a  special  form 
for  the  Poisson  parameter  \(t),  and  taking  A(u)  du  =  t-zb(t) 

Hasofer  was  able  to  write  the  emptiness  probabilities  for  the 

00 

queue  at  time  t  as  a  power  series  in  z,  P(x,t)  =  ^  znF  (t). 

n=0  n 

Using  Reich's  equations,  some  complex  analysis  and  the  added 
assumption  that  b(t)  and  b'(t)  are  uniformly  bounded,  Hasofer 
was  able  to  find  a  general  expression  for  the  Laplace  trans¬ 
form  of  the  F  . 

n 

More  concise  results  were  obtained  when  the  Poisson  pa 

n 

rameter  was  assumed  to  be  periodic.  For  A(t)  =  t-z  Y'  a 

n=l  n 

sin(nwt+<b  )  ,  Hasofer  showed  that  the  functions  F  have  the 
n  n 

asymptotic  form: 

oo 

Fn(t)  =  Y2  lAKncos  ^Kwt)  +  BKnsin(Kwt)] 

K-l 

and  the  Laplace-Stielt jes  transform  of  the  waiting  time  has  a 
similar  asymptotic  form. 

Leesc  and  Boyd,  in  a  later  work  [LB  1966],  revisited 
the  numerical  method  proposed  by  Luchak  and  showed  that  it 
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becomes  intractable  if  the  system  of  interest  is  observed  for 
any  significant  period  of  time.  They  also  analyzed  a  number 
of  other  approaches  ranging  from  direct  solutions  of  the  Kol¬ 
mogorov  difference-differential  equations  and  generating  func¬ 
tion  techniques,  to  a  Taylor  series  expansion  of  the  queue 
length  probabilities,  a  matrix  approach  using  step  function 
approximations  of  the  arrival  parameter  and  an  integral  equa¬ 
tion  technique  developed  by  Wragg  [WRAG  1963]  which  is  somewhat 
similar  to  that  described  by  Clarke.  in  all  cases,  based  on 
computer  resources  available  at  that  time,  the  computational 
requirements  were  found  to  be  excessive. 

Neuts  [NEUT  1970]  examined  the  transient  behavior  of 
two  queues  in  tandem,  where  the  second  queue  is  assumed  to  have 
finite  buffer  capacity.  The  input  to  the  first  queue  is  assumed 
to  be  Poisson,  while  the  service  times  of  the  first  queue  are 
assumed  to  be  derived  from  a  general  distribution.  The  analysis 
relies  on  several  embedded  Markov  renewal  processes. 

Rider  (RIDE  1976]  described  an  approximation  to  the  ex- 
oerteJ  queue  length  of  a  time-dependent  M/M/1  queue.  He  reduces 
'  f'o  difference-differential  equations  associated  with  the  M/M/1 
jueue  to  a  single  equation  for  the  expected  queue  length  which 
is  dependent  on  P  (t),  the  probability  that  the  system  is  empty 
a i  time  t.  Exact  expressions  for  PQ(t)  under  restricted  condi¬ 
tions  ar^  provided  as  well  as  an  approximation  relating  PQ(t) 
to  tne  expected  queue  length  at  time  t  for  those  cases  when  the 
'-ate  of  change  of  the  queue  length  is  less  than  the  service  data 

Ross  | ROSS  1978]  examined  a  non- stat ionary  M/G/l  queue 
wierc  the  t i me -dependent  arrival  parameter  A(t)  is  itself  gov- 
•'ned  by  a  stochastic  process.  Ross  conjectures  that  the  ex- 
pec ted  time  a  customer  spends  in  this  queue  is  greater  than  the 
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expected  waiting  time  for  a  customer  in  a  queue  whose  arrival 
process  parameter  is  stationary. 


Kotiah  IKOTI  1978]  proposed  approximations  for  the  time' 
dependent  expected  queue  length  of  an  M/M/1  queue.  He  discusses 
iterative  approaches  that  yield  rational  approximations  to  the 
complex  root  U*  which  lies  inside  the  unit  circle  of  a  quadratic 
obtained  by  a  Laplace  transform-generating  function  technique. 
Extension  of  the  methods  to  the  queues  M/M/2  and  M/E^/l  are  also 
discussed . 

Middleton  [M1DD  1979]  also  examined  the  numerical  solu¬ 
tion  of  inverse  transform,  and  established  a  number  of  tables 
which  can  be  used  in  the  solution  of  time-dependent  M/G/l  queues 

McClish  ]MCCL  1979}  examined  the  M(t)/G/1  queue  where 
it  was  assumed  that  the  arrival  parameter  is  of  the  form  A(t)  = 
Xq  (1  +  e$(t>),  and  the  formulas  are  series  expansions  in  e. 

It  is  shown  that  under  a  variety  of  circumstances,  periodicity 
of  the  input  is  sufficient  to  guarantee  that  a  quasi-limiting 
distribution  of  the  queue  size  exists. 

Kambo  and  Bhalaik  I  KB  1979]  examined  two  M/M/»  queues 
in  tandem  without  feedback,  where  the  arrivals  and  departures 
have  time-dependent  parameters.  They  provide  necessary  and 
sufficient  conditions  for  the  two  queue  length  distributions  to 
be  independent.  They  also  observe  that  with  certain  restric¬ 
tions  on  the  service  rate,  the  system  is  Poisson  in  the  limit¬ 
ing  case. 


Agrawala  and  Tripathi  (AT  1979]  present  a  method  of 
obtaining  the  transient  solution  of  a  general  single  server 
queue  which  relies  on  characterizing  the  queueing  system  in 
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terms  of  the  virtual  waiting  time.  The  approach  makes  no  as¬ 
sumptions  about  the  independence  of  the  arrivals  and  provides 
exact  expressions  for  the  virtual  waiting  time  distribution  and 
expected  virtual  waiting  time. 

An  example  of  the  application  of  transient  analysis  to 
a  real  problem  is  described  by  Bookbinder  and  Martell  IBM  1979]. 
They  employ  a  finite  buffer,  time-dependent  single  server  queue 
to  model  helicopter  allocation  for  forest  fires.  A  straight¬ 
forward  Runge-Kutta  method  is  used  to  solve  the  associated  finite 
set  of  difference-differential  equations. 

One  final  continuous  time  technique  of  note  involves 
the  diffusion  approximation  ([NEWE  1968]  [NEWE  1971]).  The 
technique  approximates  queue  processes  and  characteristics  by 
continuous  functions,  but  is  accurate  only  for  heavily  utilized 
(p  near  1)  queues. 

Of  the  techniques  employing  a  discrete  time  approach, 
one  of  the  first  was  offered  by  Galliher  and  Wheeler  [GW  1958]. 

In  their  technique,  a  system  is  observed  at  fixed,  equidistant 
points  in  time.  The  number  of  customers  in  the  system  at  each 
of  these  instants  is  examined  via  a  Markov  chain  analysis  whose 
transition  probabilities  are  derived  from  a  Poisson  distribution. 

Leese  and  Boyd  [LB  1966]  also  employed  a  discrete 
time  approach  to  study  the  M(t)/M/1  queue.  Their  technique 
involves  the  juxtaposition  of  a  series  of  M/M/1  queues;  a  fi¬ 
nite  time  analysis  then  has  to  be  performed  during  each  in¬ 
terval  in  which  the  arrival  parameter  is  constant.  The  solu¬ 
tion  they  proposed,  however,  involves  the  use  of  infinite  sums 
of  Bessel  functions  to  derive  the  transient  characteristics  of 
each  queue,  and  quickly  becomes  intractable. 
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Eisen  and  Tainiter  [ET  1963]  and  Yechiali  and  Naor 
[YN  1969]  investigated  queues  for  which  the  arrival  (and  pos¬ 
sibly  service)  rate  is  not  a  deterministic  function  of  time. 

The  arrival  rate  is  instead  assumed  to  be  a  heterogeneous  Pois¬ 
son  process  governed  by  an  extraneous  phase  process  which  is 
itself  a  continuous  time  Markov  chain.  Customers  arrive  in  a 
Poisson  stream  with  the  arrival  rate  dependent  upon  the  phase 
of  the  queue;  the  amount  of  time  spent  in  each  phase  is  expo¬ 
nentially  distributed.  Service  times  are  either  constant  or 
vary  with  each  phase.  The  major  results  derived  are  for  the 
case  of  a  process  with  two  phases,  where  the  service  time  is 
exponentially  distributed. 

Neuts  |NEUT  1971a]  extended  the  work  of  Yechiali  and 
Naor  to  examine  a  sequence  of  M/G/l  queues  each  having  param¬ 
eters  which  are  fixed  over  time  intervals  whose  lengths  are 
exponentially  distributed  and  where  the  parameters  attain  new 
values  in  accordance  with  a  given  set  of  probabilities.  This 
approach  is  generally  inapplicable  to  most  useful  queueing 
problems  where  significant  parameter  changes  occur  in  a  more 
deterministic  fashion. 

Another  technique  proposed  by  Neuts  [NEUT  1971b]  for 
the  GI/G/1  queue  consists  of  dividing  the  time  axis  into  units 
which  correspond  to  the  shortest  possible  service  time  and 
estimating  the  distribution  of  arrivals  during  each  such  in¬ 
terval.  A  sequence  of  queue  length  distributions  can  then  be 
derived  by  studying  a  Markov  chain  whose  states  are  2-tuples 
comprised  of  the  number  of  customers  in  the  system  and  the 
number  of  units  of  service  remaining  for  the  customer  currently 
being  served. 

More  recently,  generalized  version  of  the  methods 
initially  developed  by  Luchak  and  Clarke  have  appeared.  Koopman 
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(KOOP  1972]  studied  a  system  which  has  a  periodic  arrival  proc¬ 
ess  and  a  finite  queue.  The  arrival  and  service  time  parameters 
are,  therefore,  dependent  both  upon  time  and  the  current  queue 
length.  Koopman's  technique  then  involved  a  numerical  solution 
of  a  finite  set  of  differential  equations. 

Moore  [MOOR  1972]  (MOOR  1975]  employed  an  embedded 

X  Y 

Markov  chain  approach  to  solve  the  M  (t)/E  /I  queue,  where  the 

Y 

notation  M  denotes  a  compound  Poisson  process  for  which  cus¬ 
tomers  arrive  in  groups  of  random  size  X  which  are  exponentially 
distributed.  Since  the  compound  Poisson  process  permits  the 
ratio  of  the  variance  to  the  mean  to  be  greater  than  1,  a  wider 

range  of  input  processes  can  be  approximated.  The  generality 
v 

of  the  E  service  distribution  is  well  known  and  was  established 
in  [CAVE  1954]  [LUCH  1956].  Moore's  basic  approach  uses  the 
Chapman-Kolmogorov  equations  for  an  M/G/l  queue.  The  regenera¬ 
tion  points  for  the  imbedded  Markov  chain  are  the  departure 
instants  of  customers  from  the  queue. 

In  a  later  paper,  Minh  (MINH  1978]  explores  the  M^(t)/G/1 
queue  in  discrete  time.  In  this  technique,  the  system  is  ob¬ 
served  at  equally  spaced  intervals  and  all  events  of  the  system 
{e.g. ,  arrivals,  transfers  from  queue  to  service,  and  depar¬ 
tures)  are  assumed  to  occur  at  instants  immediately  prior  to 
these  epochs.  Minh  obtains  the  Chapman-Kolmogorov  difference 
equations  for  a  multivariate  Markov  chain  which  are  then  used 
to  express  the  vitual  waiting  time  probabilities  in  terms  of 
the  emptiness  probabilities.  A  recurrence  relation  for  calcu¬ 
lating  the  emptiness  probabilities  is  also  derived.  The  tech¬ 
nique  allows  the  calculation  of  such  measures  as  expected  de¬ 
parture  and  waiting  times  for  each  customer,  expected  number 
of  customers  in  the  system  at  each  epoch,  and  residual  service 
time  probabilities. 
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It  has  become  apparent  that  the  most  viable  techniques 
for  analyzing  the  transient  behavior  of  queues  center  about  a 
discrete  time  approach  ((MOOR  1975]  [MINH  1978]).  The  rigid 
assumptions,  mathematical  intractability  and  computational 
complexity  inherent  to  continuous  time  techniques  limit  their 
flexibility  and  extensibility.  Yet,  while  the  discrete  time 
approaches  reduce  (if  not  overcome)  these  difficulties,  they 
still  present  sizable  challenges  in  trying  to  accurately  rep¬ 
resent  queues  via  discretely  time-varying  processes  and  ana¬ 
lyzing  the  behavior  of  complex  networks  of  such  queues;  in 
fact,  little  work  has  been  performed  to  date  in  the  latter 
area . 


Our  objective  in  this  paper,  then,  is  to  begin  to  lay 
the  groundwork  for  the  approximate  (yet  effective)  transient 
analysis  of  queueing  networks.  We  propose  to  do  this  in  two 
stages.  In  Section  2.1,  we  will  develop  a  modified  version  of 
the  finite  time  analysis  initially  proposed  by  Leese  and  Boyd 
(LB  1966]  for  the  M(t)/M/1  queue.  )ur  modification  centers 
about  providing  an  approximate  method  (based  on  work  by  Stern 
(STER  1979])  for  deriving  the  transient  behavior  of  an  M/M/1 
queue  which  will  replace  the  exact  solution  technique  for  ob¬ 
taining  the  transient  behavior  based  upon  infinite  sums  of 
Bessel  functions.  This  technique  is  both  computationally  ef¬ 
ficient  and  can  yield  results  to  any  desired  level  of  accuracy. 

In  Section  2.2  we  will  show  that  by  making  certain 
assumptions  about  the  ergodicity  of  the  input  and  service  proc¬ 
esses  of  the  M(t)/M/1  queue  over  finite  intervals,  an  input 
process  which  is  discretely  time-varying  exponential  ti.e., 

M(t))  yields  a  departure  process  which  is  also  of  the  form 
M(t).  This  result,  which  is  an  approximation,  greatly  facili¬ 
tates  the  integration  of  a  set  of  M(t)/M/1  queues  into  a  network, 
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and  the  subsequent  transient  solution  of  that  network.  It  is 
worth  noting  that  to  date,  only  one  result  relating  to  depar¬ 
ture  processes  of  time  dependent  queues  has  been  derived;  this 
when  Mirasol  [MIRA  1963]  and  Kendall  [KEND  1964]  showed  that 
the  departure  process  for  the  M(t)/G/»  queue  is  M(t). 

Following  these  two  developments,  the  results  of  sev¬ 
eral  experiments  validating  our  analysis  will  be  provided  in 
Section  3.  A  summary  of  our  current  work  and  a  description  of 
future  research  are  developed  in  Section  4. 


2.  APPROXIMATE  TRANSIENT  ANALYSIS  OF  THE 

M(t)/M/1  QUEUE 

2.1  THE  APPROXIMATION  TECHNIQUE 

The  precise  situation  with  which  we  are  concerned  is 
shown  in  Fig.  2-1.  In  consonance  with  other  discrete  time  ap¬ 
proaches,  we  will  make  the  assumption  that  the  arrival  rate 
parameter,  A(t),  varies  (or  can  be  approximated  as  varying) 
with  time  in  a  step-wise  fashion,  and  that  between  those  in¬ 
stants  when  A ( t )  changes,  it  can  be  assumed  to  be  stationary. 
We  will  further  assume  that  the  service  parameter  remains 
constant  for  all  intervals. 


A^  »M  A2  »M  *3  »H 

| - I - 1 - - 

t0  tl  l2 

Figure  2-1 

The  instants  at  which  A(t)  changes  are  given  by  tp, 
t-p  ...  and  the  intervals  bounded  by  those  instants  are  1^  = 
(t(),tl>,  *2  =  etc*  Since  the  service  parameter  remains 

constant  for  all  intervals,  we  have  within  each  interval  an 
M/M/1  queue.  Our  analysis,  therefore,  conveniently  decomposes 
into  the  sequential  analysis  of  a  set  of  (stationary)  M/M/1 
queues  over  finite  intervals;  it  is  the  analysis  of  these 
individual  queues  which  concerns  us  in  this  section. 
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An  approximate  transient  analysis  of  the  stationary  M/M/1 
queue  has  been  performed  by  Stern  (STER  1979]  and  relies  heavily 
upon  techniques  developed  by  Keilson  [KEIL  1974]  [KEIL  1964]  (KEIL 
1965]  [KEIL  1966].  In  trying  to  overcome  the  complexity  associ¬ 
ated  with  deriving  exact  expressions  for  such  quantities  as  the 
transient  expected  queue  length  or  throughput,  which  involve  in¬ 
finite  sums  of  Bessel  functions  [SAAT  1961],  Stem  approximates 
the  M/M/1  queue  by  a  finite  buffer,  M/M/l/K  queue.  The  justifi¬ 
cations  for  such  an  approximation  include  the  facts  that: 


•  Over  a  finite  interval,  only  a  limited  num¬ 
ber  of  arrivals  can  occur;  hence,  providing 
infinite  buffer  capacity  is  unnecessary. 

•  Few  real  systems  have  infinite  buffer  space 
and  consequently  do  not  exhibit  true  M/M/1 
behavior . 


Stern  begins  with  the  fundamental  transition  rate 

equation 


p'  =  Q  *  P  (1) 

where  Q  is  the  (K+1)*(K+1)  infinitesimal  transition  rate  matrix 
associated  with  the  finite  Markov  chain  underlying  the  M/M/l/K 
queue,  p  is  the  state  probability  vector,  and  p'  is  the  deriva¬ 
tive  of  p.  By  changing  variables,  it  is  possible  to  obtain 

u'  =  E**(-l/2)  *  p  (2) 

where 
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and 


en  =  p**n  (1-P)/<1-P**(K+1>> »  P  =  A/M 


which  is  just  the  steady  state  solution  of  being  in  state  n 
(<K)  for  the  M/M/l/K  queue. 


Proceeding,  we  can  define 


where 


u '  =  -p  *  A  *  u 


p  Vp  0  ... 

-Vp  (1+p)  -  Jp  0  ... 


~ylp  1 


(3) 


Now,  the  transformation  in  (2)  symmetrizes  the  matrix  Q  in 
(1).  Consequently,  the  finite  Markov  chain  associated  with  Q 
is  time  reversible  and  a  spectral  representation  of  the  time 
dependent  behavior  of  the  chain  is  available  for  which  the 
eigenvalues  and  eigenvectors  are  real.  In  fact,  both  Q  and  A 
have  real  eigenvalues  and  the  solution  to  equation  (3)  may, 
therefore,  be  written  as 


u(t)  =  exp(-p*l.*t)  u.  uT  u(0)  (4) 

Li=0  1  11  J 

where  1  and  ui  are  the  i^  eigenvalue  and  eigenvector  (ar¬ 
ranged  in  ascending  order)  of  A,  respectively.  It  is  worth 
noting  that  1  =0  and  l.>0,  i*0 . 

O  1 
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The  transient  probability  distribution  for  the  number 
in  system  can  then  be  given  by 


p(t)  =  p 


IV 


c  .  p . 

1*1 


exp( -p*l^*t) 


(5) 


where  pQ  is  the  steady  state  queue  length  distribution,  and 

Pi  =  E**(l/2)  *  u.  (5.1) 

Ci  =  u?  *  E**(-l/2)  *  p( 0 )  (5-2) 


The  expected  queue  length  at  time  t  is  seen  to  be  given  by 


n(t)  =  pT*p(t)  ,  p  =  [0,1,2 . K]T 

or 


K  T 

n(t)  =  nQ  +  ci  p  pt  exp(-p*li*t) 


(6) 


where  nQ  is  the  steady  state  expected  queue  length. 

The  application  of  (5)  and  (6)  to  our  particular  prob¬ 
lem  is  immediate.  Since  we  have  assumed  stationarity  within 
each  of  the  intervals  I ^ ,  we  are  dealing  with  a  sequence  of 
M/M/1  queues  which  we  can  approximate  as  a  sequence  of  M/M/l/K 
queues.  Some  fundamental  questions  remain,  however.  The  first 
concerns  the  way  in  which  the  behavior  of  the  queue  is  approxi¬ 
mated  when  p  >  1.  The  technique  developed  by  Stern  is  valid 
only  when  p  <  1,  since  the  associated  Markov  chain  is  otherwise 
not  ergodic  and  the  time  reversibility  argument  does  not  apply. 
As  a  result,  we  will  employ  the  following  alternative  approxima 
tion.  If  we  assume  that  the  queue  length  at  the  beginning  of 
some  interval  K  for  which  p  >  1  is  given  by  n(tj_^),  then  we 
will  approximate  the  transient  expected  queue  length  by: 
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n(  t  )  =  n(tj_j)  +  (A  -  p)  Mt-tj^)  ,  I  e  Ij 

A  second  question  centers  about  the  selection  of  the 
value  of  Kj  in  a  particular  interval  Ij  (for  which  p  <  1). 

Our  objectives  in  selecting  a  value  for  Kj  are  threefold: 


•  To  ensure  that  the  limiting  (i.e.,  steady 
state)  expected  queue  lengths  of  our 
M/M/l/K  approximation  and  the  ideal  M/M/1 
queue  are  within  a  sufficient  distance 

of  each  other. 

•  To  ensure  that  the  rate  at  which  our 
M/M/l/K  approximation  approaches  its 
steady  state  value  is  comparable  to  that 
associated  with  the  M/M/1/®  queue.  The 
term  commonly  used  to  quantify  this  rate 
of  approach  to  steady  state  is  the  relaxa¬ 
tion  time. 

•  To  ensure  that  computational  tractability 
is  maintained. 


In  this  regard,  two  alternative  methods  of  selecting  Kj  were 
considered.  The  first  sought  to  define  Kj  based  on  a  2o  esti 
mate  of  the  queue  length  during  the  interval  I j .  Thus,  if  at 
the  beginning  of  Ij ,  n{tj_^)  =  0,  then  Kj  would  be  defined  as 

K .  =  n  .  +  2o 
J  oj  n 

J  J  oj 

where 

n  .  =  p ./( 1 -p  .  ) 
oj  pj/v  pj' 


the  steady  state  expected  queue  length 


and 


2o  -  2  /pW)/(l-p,)3 


oj 
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=  2p  j/ (1-p  j )  V(l/Pj  +  Pj>/^  -  Pj^ 

=  twice  the  standard  deviation  of  n 

oj 


On  the  other  hand,  if  at  the  beginning  of  an  interval 
I j ,  n(tj_^)  >  0,  a  choice  for  Kj  of 


Kj  =  "‘/j-/  +  "oj  +  2an 


OJ 


will  do  the  same. 


A  second  method  of  defining  Kj  is  closely  tied  to  the 
relaxation  time  of  the  associated  approximation.  By  Stern’s 
method,  the  relaxation  time  of  a  particular  M/M/l/K  approxima¬ 
tion  is  given  by: 


t  s  1/p  (1+p  -  2,/p  cos  0) 


where 


and  the 


0  =  n/K+1 
expression 


1+p  -  2 Vp  cos  0 


is  equal  to  the  least  non-zero  eigenvalue  associated  with  the 
matrix  A  defined  earlier.  Obviously,  as  K  •*  ®,  1  approaches 
its  minimum  value  (and  x  its  maximum). 


It  was  decided  that  Kj  would  be  defined  such  that  the 
relaxation  time  associated  with  the  M/M/l/K  queue  was  equal  to 
2/3  that  of  the  M/M/1  case.  Thus,  if  at  the  beginning  of  I j , 
n ( t j  _  j  )  =  0,  then 
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,1  (  M  4  ) 

K.  =  INT  (cos  i(l . 5-  .25  - J-))  +1 

J  *Aj 


where  INT  is  a  function  which  yields  the  greatest  integer  less 
than  or  equal  to  its  argument.  If  at  the  beginning  of  I.,  n 
(tj _^)  >0,  then  we  define 


Kj  =  INT  (n  (tj_1)  +  cos 


-1 


(1.5-  .25 


(M+X ,  ) 


J-  ))  +1  (7) 


Tables  2-1  and  2-2  show  how  the  M/M/l/K  approximation 
associated  with  the  first  and  second  definitions,  respectively, 
of  Kj  compare  with  the  M/M/1  case  and  each  other.  The  initial 
queue  length  is  assumed  to  be  zero  in  all  cases.  As  can  be 
seen,  it  is  not  until  p>.88  that  the  2o-based  definition  of 
yields  values  for  the  expected  steady  state  queue  length  and 
relaxation  time  which  are  closer  to  the  exact  M/M/1  case  than 
those  of  equation  (7).  As  a  consequence,  equation  (7)  was 
chosen  to  serve  as  our  general  definition  for  . 


The  third  question  stems  from  the  fact  that  the  value 
of  Kj  usually  varies  over  different  intervals,  creating  a  prob¬ 
lem  for  the  effective  approximate  analysis  of  a  series  of  such 
queues.  From  (5.2)  and  (7),  it  is  evident  that  the  distribu¬ 
tion  p(t)  plays  an  important  part  in  the  transition  between 
distinct  intervals.  In  fact,  the  initial  queue  length  distri¬ 
bution  p(0)  appearing  in  equation  (5.2)  for  a  particular  in¬ 
terval  will  be  set  equal  to  the  distribution  p(t)  at  the  end 
of  the  preceeding  interval.  It  will  therefore  be  necessary  to 
translate"  the  + ^ -d imensional  queue  length  distribution 
vector  p j ( t  )  a*  the  end  of  some  interval  1^  into  the  K^  +  ^  +  l- 
dimensional  queue  length  distribution  vector  Pj+j^O)  at  the 
beginning  cf  the  succeeding  interval  Ij  +  j. 
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TABLE  2-1 


2a -BASED  DEFINITION  RESULTS 


p 

K. 

J 

_  -k 

N 

APP 

- 

N 

*  EXACT 

napp/nexact 

0.10 

0 

0.11111 

0.20 

1 

0.16667 

0.25000 

0.66667 

0.30 

2 

0.34532 

0.42857 

0.80576 

0.40 

3 

0.56158 

0.66667 

0.84236 

0.50 

5 

0.90476 

1.00000 

0.90476 

0.60 

8 

1.40838 

1.50000 

0.93892 

0.70 

14 

2.26178 

2.33333 

0.96933 

0.80 

29 

3.96282 

4.00000 

0.99070 

0.81 

32 

4.23161 

4.26316 

0.99260 

0.82 

35 

4.52711 

4.55556 

0.99376 

0.83 

38 

4.85510 

4.88235 

0.99442 

0.84 

42 

5.22614 

5.25000 

0.99545 

0.85 

47 

5.64701 

5.66667 

0.99653 

0.86 

52 

6.12496 

6.14286 

0.99709 

0.87 

59 

6.67820 

6.69231 

0.99789 

0.88 

67 

7.32192 

7.33333 

0.99844 

0.89 

77 

8.08211 

8.09091 

0.99891 

0.90 

89 

8.99314 

9.00000 

0.99924 

0.91 

105 

10.10628 

10.11111 

0.99952 

0.92 

126 

11.49680 

11.50000 

0.99972 

0.93 

155 

13.28382 

13.28571 

0.99986 

0.94 

196 

15.66567 

15.66667 

0.99994 

0.95 

259 

18.99958 

19.00000 

0.99998 

0.96 

363 

23.99986 

24.00000 

0.99999 

0.97 

560 

32.33331 

32.33333 

1.00000 

0.98 

1029 

49.00000 

49.00000 

1.00000 

0.99 

2899 

99.00001 

99.00001 

1.00000 

^Expected  steady  state  queue  length  for  this  ap¬ 
proximation  . 

**Expected  steady  state  queue  length  for  M/M/1 
queue . 


I 
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TABLE  2-2 


[ 

P 

! 


RELAXATION  TIME-BASED  DEFINITION  RESULTS 


K 


N 


APP 


N 


EXACT 


^APP7^  EXACT 


0.10 

0.20 

0.30 

0.40 

0.50 

0.60 

0.70 


3 

5 

7 

9 

12 

17 

24 


0.11071 

0.24962 

0.42805 

0.66562 

0.99841 

1.49817 

2.32998 


0.11111 

0.25000 

0.42857 

0.66667 

1.00000 

1.50000 

2.33333 


0.99640 

0.99846 

0.99878 

0.99843 

0.99841 

0.99878 

0.99856 


0.80 

0.81 

0.82 

0.83 

0.84 

0.85 

0.86 

0.87 

0.88 

0.89 


39 

42 

44 

47 

50 

54 

58 

63 

69 

76 


3.99468 

4.25816 

4.54960 

4.87609 

5.24299 

5.65945 

6.13480 

6.68360 

7.32424 

8.08115 


4.00000 

4.26316 

4.55556 

4.88235 

5.25000 

5.66667 

6.14286 

6.69231 

7.33333 

8.09091 


0.99867 

0.99883 

0.99860 

0.99872 

0.99866 

0.99873 

0.99860 

0.99871 

0.99876 

0.99879 


0.90 

0.91 

0.92 

0.93 

0.94 

0.95 

0.96 

0.97 

0.98 

0.99 


84 

94 

106 

122 

143 

173 

217 

291 

439 

884 


8.98902 

10.09890 

11.48572 

13.26937 

15.64722 

1897685 

23.97024 

32.29327 

48.93933 

98.87863 


9.00000 

10.11111 

11.50000 

13.28571 

15.66667 

19.00000 

24.00000 

32.33333 

49.00000 

99.00001 


0.99878 

0.99879 

0.99876 

0.99877 

0.99876 

0.99878 

0.99876 

0.99876 

0.99876 

0.99877 


A  problem  arises,  however,  when  for  two  contiguous 
intervals,  Ij  and  Ij  +  ^,  ^j^j+1"  *n  t^1^s  case »  the  dimensions 
of  the  associated  Pj(t)  and  p^+^(t)  are  different,  and  an  exact 
assignment  cannot  be  made.  Some  possible  approaches  for  solving 
this  problem  include: 


1)  Setting  all  K.'s  equal  to  a  number  which  is  greater 
than  or  equalJthe  largest  possible  value  of  any 
K.  over  the  intervals  of  interest.  The  advantage 
oi  this  approach  is  that  it  guarantees  that  all 
p.(t)'s  have  the  same  dimension.  The  disadvan¬ 
tages,  however,  include: 


•  Computational  inefficiency,  since  one  is 
always  working  with  the  maximum  M/M/l/K 
approximation  between  intervals. 


•  It  will  usually  be  impossible  to  deter¬ 
mine  just  what  the  appropriate  maximum 
value  for  K  should  be.  Recall  that  the 
K's  are  determined  dynamically  as  the 
analysis  proceeds  and  rely  on  information 
obtained  not  only  from  the  current  inter¬ 
val  but  the  preceding  interval  as  well. 


2)  If  Kj  £  Kj+1’  Pr°jecting  (in  the  strict  mathematical 
sense)  p^Ct^-t^ _^)  directly  in  Pj+^(0).  For  example, 

Pj,l(t) . Pj,m(t))  and 

p..,  „(t))  where  m<n,  then 
j+l>n  - 

Pj+1(t)  implies  that  p^  Q(t)  = 


if  we  let  Pj(t)  =  (Pj  q( t ) » 

pj+1<t>  =  (Pj*l>o<t>’- 
projecting  p^(t)  into 


pj*l,0U)'  pj,l(t>  =  Pjtl,l(t)’  ••••  pj,m(t)  *  pj+l,k(t) 
and  Pj+^,£(t)  =  0  V£  such  that  k<£<n.  On  the  other 

hand,  if  K.>K..,,  then  we  must  "scale"  the  dimension 

of  Pj(t)  down  to  that  of  Pj+^(t)  ensuring,  however, 

that  n(tj)  =  E (Pj ( t j - 1 j _^ ) ]  equals  n(tj)  =  E[pj+1(0)]. 


One  possible  method  of  performing  this  scaling  is  as 

follows.  Let  nmax  be  the  least  integer  which  is  greater  than  or 

equal  to  n(t.).  Define  p^  +  1  Q(0)  =  1  -  n(tj)/nmax,  Pj+i>nmax<°) 

=  n(t.)/nmax,  and  p.  ,  .  (0)  =  0  for  all  i/0  or  nmax.  It  is 

J  J  A  > 
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obvious  then  that  with  this  representation, 


KJ+1 


Z_Q  pj+l,i<0>  =  1  and  ElPj+i(°>)  = 


and  we  have  achieved  the  required  mapping  between  p^.  and  Pj  +  ^ 


Using  Stern's  basic  technique  together  with  the  solu¬ 
tions  described  above,  a  useful  and  computationally  efficient 
method  becomes  available  for  analyzing  the  approximate  tran¬ 
sient  behavior  of  both  a  stationary  M/M/1  queue  and  the  (dis¬ 
cretely  time- varying)  M(t)/M/1  queue.  In  the  next  section,  we 
investigate  those  aspects  of  the  M(t)/M/1  queue  which  will 
permit  us  to  extend  the  analysis  described  above  to  networks 
of  such  queues. 


2.2  APPROXIMATE  CHARACTERIZATION  OF  THE  DEPARTURE  PROCESS 

This  section  develops  an  approximate  characterization 
of  the  departure  process  for  the  M(t)/M/1  queue  described  above 
and  details  how  this  characterization  can  potentially  facilitate 
the  solution  of  general  networks  comprised  of  such  queues.  Al¬ 
though  several  papers  [CD  1974]  [NATV  1975]  [DALE  1976]  have  ex¬ 
amined  the  topic  of  departure  processes  for  a  wide  variety  of 
queues,  none  have  discussed  in  any  depth  the  nature  of  the  de¬ 
parture  process  during  a  transient  period. 

In  order  for  us  to  begin  to  characterize  the  departure 
process,  we  make  the  following  observations  and  approximations: 

•  Using  the  techniques  described  in  Section 
2.1,_we  can  obtain  an  approximate  value 
for  n(t)  for  any  t.  Given  that  each 
customer  has  a  mean  service  time  of  1/p , 
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the  expected  virtual  waiting  time  for 
the  queue  at_the  beginning  of  some  inter¬ 
val  I..-,  is  n(t.)/p.  This  means  that 

j+1  j 

during  the  interval  [t^,  t j+(n( t ^ )/p ) )  = 

[tj,tj+v>,  the  server  is  expected  to  be 

busy.  The  departure  process  will,  there¬ 
fore,  appear  as  the  service  process  and 
have  parameter  p. 

•  Given  that  the  expected  virtual  waiting 
time  at  the  beginning  of  some  interval 
Ij+2  is  non-zero  and  that  Aj+^/0,  then 

the  departure  process  will. have  param¬ 
eter  p  up  to  the  instant  T1  when  the 
server  first  goes  idle  (which  may  exceed 
t j+v  above ) . 

•  After  the  first  instant  T1 ,  t^+v  i.  Tx  < 

tj+^,  such  that  the  server  goes  idle,  we 

approximate  the  departure  process  as  the 

input  process  with  parameter  This 

approximation  makes  liberal  use  of  Burke's 
Theorem  [BURK  1956]  which  states  that  M 
=>  M;  its  accuracy  will  depend  in  part 
upon  the  relaxation  time. 


Our  first  objective  then  is  to  determine  when  (and 
if)  the  first  (expected)  instant,  Tj ,  occurs  in  interval  Ij 
such  that  the  server  goes  idle.  We  will  assume  that  Pj(tj_^) 
and  n(tj^)  are  known.  While  these  last  parameters  are  as¬ 
sumed  to  be  derived  by  the  methods  presented  in  Section  2.1, 
our  subsequent  analysis  deals  explicitly  with  the  M/M/1  queue 
and  does  not  involve  its  M/M/l/K  approximation. 


Recalling  that  fundamental  to  any  M/M/1  queue  is  a 
birth-death  process  with  transition  rates  A  and  p,  we  find 
ourselves  dealing  with  an  infinite  Markov  chain  whose  states 
are  just  the  number  in  system.  1 f  we  denote  the  passage  time 
probability  density  function  (pdf)  from  a  state  n  to  state  0 
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by  sn  q ( t ) ,  and  the  corresponding  random  variable  by  Tn  q,  it 
is  immediate  that 

s  =  sn(t)  *  Sn-1<T>  *  •••  •  sj(t)  (8) 

n ,  0 

where  s-(t)  =  ,  (t)  and  6  is  the  convolution  operator.  In 

n  n ,n-l 

terms  of  random  variables  we  have 

t  =  Tr,  +  T„  l  „  9  +  •  •  •  +  T1  n 

rn  q  n ,n  1  n* 1 ^n*2  1,0 

A  recursive  probabilistic  argument  obtains  s  aS 

follows.  Let  v  =  A+p  .  If  we  assume  that  we  are  in  sBihe  state 

n,  the  dwell  time  in  that  state  has  a  pdf  given  by  v*exp(-v*t). 

Furthermore,  with  probability  A/A+p  there  will  be  a  transition 

to  state  n+1 ,  while  with  probability  p/A+p  the  transition  will 

be  to  state  n-1.  Hence, 


s-<t) 


V  ve" 


VT  +  £  ve‘VT  •  Sn+1(T)  •  sn(T) 


The  recursion  is  obtained,  apart  from  inversion  difficulties, 
by  using  Laplace  transforms.  The  transform  of  (9)  is  then 
seen  to  be 


°n<s> 


i+Vff  +  S5X+?  '  %+l<s)  ’ 


s+A+p 


11  -  '  °S+l(s>> 


s+A+p -  Ao-+^(s' 


from  which  one  can  further  derive 


d  ,  .  -l,a  ■  x  dl  t’S+l(s)> 

-T—  O  -  (  S  )  - - A  - 

4  kl  n  (s+A+p  -  Aa-  +  1(s) 
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Clearly , 


-  o- (0)  =  E [T- ]  =  T- 
ds  n'  1  n  n 


(12) 


where 


T-  =  T  , 
n  n,n-l 


Hence , 


T-  =  m_1(1  +  *  T-+1 ) 


(13) 


The  mean  time  from  state  n  to  n-1  may  be  described  in 
more  compact  form  with  the  aid  of  the  parameters  nR  defined 


nQ  =  1  ;  nn=  U/p)n  =  pn 


where 


Arc  =  un 
n  M  n+1 


Multiplying  (13)  *>y  p/r 


,  we  obtain 


p 7i  T-  -  p 7i  ,,1-  ,  =  n 
n  n  n+1  n+1  n 


from  which  we  can  derive 


T-  =  n  . 

n  un  j 

p  n  j=n  J 


(15) 


follows  that 
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n  ,o 


1  r> 

=  m?0  j?m  "j 
n-1  ,  ®  . 

m=0  pp  J  =m 

n-1  ,  /  ,  m-1  A 

*  ?0  -is  Mp  -  £,  pJ 

m=0  MP  \  J-0  / 


n-1  , 

=  E„  A  ■  i 

m=0  pp 
n-1  , 

=  E  — 


At'a  P-X 
m=0  H 

n 

M-X 


(16) 


Thus  within  a  particular  interval  1 ^ ,  we  see  that  the  time  of 
first  idle  is  given  by: 


=  tj_2  +  n( t j  _2 )/(m-X ) 


=  t . 


j-1 


+  T_ 


n  ,o 


(17) 


Now  if  T*  >  t.  or  if  \  ./p  >  1  (i.e.,  p.  >  1),  then  it 
J  —  J  j  —  ’  j  — 

is  expected  that  the  server  will  not  be  idle  during  I j .  We 

can  therefore  make  the  following  observations: 


1)  In  the  case  where  either  >  t .  or  A  .  >  p  . 

J  -  J  J  - 

the  departure  process  during  1^  can  be 

expected  to  appear  as  the  service  process 
with  parameter  p. 
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2)  In  the  case  where  Tj  =  tj_^  (i.e.,  where 

=  0)  and  X ^  <  p ,  the  departure 

process  can  be  approximated  as  the  input 
process  with  parameter  X  ^  . 

3)  When  t._,  <  T*  <  t.  and  X.  <  p,  the  de- 

J  J  *3  J 

parture  process  can  be  expected  to  appear 
as  the  service  process  with  parameter  p 

up  to  time  after  which  is  can  be  ap¬ 
proximated  as  the  input  process  with 
parameter  X ^ . 


Based  on  these  observations,  a  general  expression  for 
the  interdeparture  time  pdf,  d(s),  within  an  interval  1^  becomes 


-p(s-t.  ,)  -p(s-t.  ,) 

d(s)  =  pj (s )pe  3  +  ( 1 -p J ( s ) )  [p jpe  J 


+  (l-pj)Xje  J 


-X.(s-t.  ,)  -p(s-t.  ,) 

J  J  1  6  pe  J  1 


] 


(18) 


where  the  density  function  pj(s)  is  defined  as: 


P  i  (  S  )  =  'j_1  ~ 

J  0  otherwise 


This  expression  really  implies,  however,  that  d(s)  is  a  step¬ 
wise  time-varying  exponential  distribution;  hence,  the  depar¬ 
ture  process  is  of  the  form  M(t).  It  is  worth  noting,  though, 
that  the  form  of  M(t)  for  the  departure  process  will  usually 
not  be  equal  to  the  form  of  M(t)  for  the  input  process.  More 
specifically,  the  instants  at  which  the  parameters  associated 
with  the  arrival  and  departure  process  change,  as  well  as  the 
values  they  acquire  between  those  instants,  will  usually  be 
di f  ferent . 
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There  are  several  implications  of  these  results. 

Probably  the  most  significant  is  the  fact  that  the  approximate 
analysis  of  networks  of  M(t)/M/1  queues  now  becomes  tractable. 
Since  the  output  of  any  M(t)/M/1  queue  is  now  in  the  same  form 
as  the  input  to  all  others,  the  interconnection  of  such  queues 
is  straightforward. 

A  second  major  implication  is  that  it  is  now  possible 
to  approximately  determine  (for  networks  of  such  queues)  just 
how  long  it  takes  a  change  in  the  arrival  rate  of  one  M(t)/M/1 
queue  to  propagate  through  a  series  of  such  queues,  and  moreover, 
when  equilibrium  will  be  re-established  for  the  entire  network. 

Finally,  from  the  point  of  view  of  potential  practical 
applications,  the  ability  to  approximately  analyze  the  transient 
behavior  of  networks  of  such  queues  lends  itself  immediately 
to  such  problems  as  the  modeling  of  dynamic  routing  strategies, 
internetworking  problems  and  network  stability  analysis. 
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EXPERIMENTAL  RESULTS 


This  section  summarizes  the  results  of  several  experi¬ 
ments  that  were  performed  to  validate  the  approximation  tech¬ 
niques  presented  in  Section  2.  In  all  cases,  the  latter  tech¬ 
niques  were  compared  to  the  output  of  simulations  or  analytic¬ 
ally-derived  exact  values.  The  experiments  had  four  objectives 


1)  Given  stationary  values  for  A  and  p  ,  to 
determine  the  relative  accuracy  of  the 
M/M/l/K  approximation  developed  in  section 
2.1  compared  to  the  exact  (analytical)  M/M/1 
case  and  simulation. 

2)  To  evaluate  the  proposed  (approximate)  method 
of  analyzing  the  behavior  of  a  noristationary 
M(t)/M/1  queue  as  a  sequence  of  stationary 
M/M/l/K  queues. 

3)  To  determine  the  validity  of  approximately 
characterizing  the  departure  process  of  an 
M ( t ) /M/1  queue  as  a  discretely  time-varying 
M(t)  process. 

A)  To  provide  a  comparison  of  the  M/M/l/K  ap¬ 
proximation  technique  with  Moore's  [MOOR 
1975]  embedded  Markov  chain  approach. 


In  sections  3.1  -  3. A  below,  each  of  these  objectives  is  sepa¬ 
rately  addressed. 


3.1  M/M/l/K  APPROXIMATION  RESULTS 

A  number  of  experiments  were  run  to  determine  the 
sensitivity  and  accuracy  of  the  M/M/l/K  approximation  described 


] 


in  section  2.1  to  variations  in  the  values  of  K,  initial  queue 
length  and  p.  Figure  3-1  illustrates  the  way  in  which  the 
transient  expected  queue  length  varies  over  time  for  K=10,  15, 

20  and  40.  In  all  cases,  the  initial  queue  length  is  assumed 
to  be  equal  to  10  and  p=.8  (with  p=50). 

The  variations  that  can  be  seen  in  the  curves  associated 
with  the  different  values  of  K  derive  from  two  sources.  The 
first  centers  about  the  fact  that  the  steady  state  values  for 
the  expected  queue  length  are  monotonically  increasing  with  K; 
hence,  the  M/M/1  value  is  asymptotically  approached  as  K  +  ». 

The  second  involves  the  relaxation  times  associated  with  par¬ 
ticular  M/M/l/K  approximations  and  the  fact  that  they  too  are 
monotonically  increasing  with  K. 

The  combined  effect  of  these  two  conditions  is  that, 
in  comparison  to  the  exact  M/M/1  case,  an  M/M/l/K  approxima¬ 
tion  based  on  a  small  value  of  K  approaches  steady  state  more 
rapidly  and  with  an  asymptotic  expected  queue  length  that  is 
less  than  an  approximation  based  on  a  larger  value  of  K.  It 
was  for  these  reasons  that,  in  attempting  to  select  a  value 
for  K  (given  a  fixed  initial  queue  length  and  p),  we  required 
that  the  M/M/l/K  approximation  provide  an  expected  steady  state 
queue  length  within  .2%  of  the  expected  M/M/1  value  together 
with  a  relaxation  time  within  33%  of  the  M/M/1  value  while  at 
the  same  time  maintaining  computational  tractability . 

Figure  3-2  illustrates  how  our  technique  performs  in 
comparison  to  simulations  based  on  the  same  underlying  para¬ 
meters  of  p=.8,  p=50  and  an  initial  queue  length  of  10.  Given 
these  conditions,  our  technique  specifies  that  we  use  an  M/M/l/K 
approximation  for  which  K=49.  As  can  be  seen,  close  agreement 
is  obtained;  in  fact,  the  average  difference  between  the  results 
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associated  with  the  approximation  and  those  associated  with 
the  simulation  was  only  0.179,  with  a  standard  deviation  of 
0.153.  Similar  agreement  was  obtained  in  other  experiments. 


3.2  M(t)/M/1  APPROXIMATION  RESULTS 

A  series  of  experiments  was  run  to  compare  the  tran¬ 
sient  expected  queue  lengths  associated  with  an  M(t)/M/1  queue 
derived  both  by  simulation  and  the  approximate  sequential 
M/M/l/K  technique  described  in  section  2.1.  Figure  3-3  illus¬ 
trates  the  outcome  of  analyzing  one  particular  M(t)/M/1  queue 
for  which  the  arrival  process  was  permitted  to  vary  such  that 
\  =  35 ,  40  and  25  for  three  consecutive  intervals,  each  of  2 
second  duration.  For  all  intervals,  pi  =50  and  was  fixed. 


Table  3-1  compares  the  expected  queue  lengths  at  the 
end  of  each  interval  obtained  by  our  approximation  and  the 
simulation  with  those  values  associated  with  an  M/M/1  queue. 
Since  for  each  interval,  the  interval  length  is  greater  than 
its  associated  relaxation  time,  the  (absolute)  difference  be¬ 


tween  the  expected  queue  length  obtained  by  simulation 

and  the  M/M/1  value  (N  ),  and  the  difference  between  N  .  and 

ss  sim 


the  expected  queue 
serve  as  relative 


length  obtained  by  our  approximation  (N 

app 

figures  of  merit. 


), 


TABLE  3-1 


END  OF  . . . 

O 

N  . 
sim 

N 

app 

N 

ss 

|N  .  -N  | 

'  sim  ss' 

|N  .  -N  | 

1  sim  app1 

j  INTERVAL 

1 

0.7 

2.308 

2.319 

2.333 

0.025 

0.011 

1 NTERVAL 

2 

0.8 

3.545 

3.818 

4.000 

0.455 

0.273 

INTERVAL 

3 

0.5 

1.018 

0.999 

1.000 

0.018 

0.019 

F igure 


To  further  demonstrate  the  "goodness  of  fit"  of  the 
approximation  results  with  those  obtained  by  simulation,  the 
table  contained  in  Figure  3-3  shows  how  the  mean  queue  lengths 
of  each  differed  over  the  three  intervals  and  includes  standard 
deviations  for  each.  Once  again,  the  close  agreement  observed 
was  repeated  in  other  similar  experiments. 


3.3  DEPARTURE  PROCESS  APPROXIMATION  RESULTS 

Several  experiments  were  run  to  validate  our  conten¬ 
tion  that  during  an  interval  whose  initial  queue  length  is 
greater  than  zero,  the  departure  process  of  an  M/M/1  queue  can 
be  approximated  as  the  discretly  time-varying  M(t)  process 
described  in  section  2.2.  Specifically,  our  contention  is 
that  between  the  beginning  of  the  interval  and  the  time  of 
first  idle,  T1 ,  the  departure  process  has  a  probability  den¬ 
sity  function  which  can  be  approximated  as  pe"^Jt  (0£t<T1). 

After  T1 ,  we  assume  that  the  departure  process  has  a  pdf  which 
can  be  approximated  by  Xe”**'  (T^t). 

Table  3-2  illustrates  how  the  expected  time  of  first 

idle  obtained  via  simulation  (T1 .  )  compares  with  our  analytic 
•  s  im 

result  (T*  )  for  p  =  .5,  .8  and  .9,  p  =  50  (and  fixed),  and 
an 

initial  queues  lengths  of  5  and  10  customers.  In  all  cases, 
the  close  agreement  predicted  by  Eqs.  16  and  17  was  observed. 

Table  3-3  shows  how  the  interdeparture  distributions 
derived  from  simulations  associated  with  several  intervals 
having  different  initial  queue  lengths  compare  to  the  predic¬ 
ted  approximate  values.  Both  the  interdeparture  distribution 
mean  and  standard  deviation  for  the  "pre-T1"  period  and  "post- 
T1"  period  are  given.  Once  again,  the  close  agreement  observed, 
combined  with  similar  results  from  other  experiments,  lend 
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TABLE  3-2 


INITIAL 
QUEUE  LENGTH 

P  (M=50) 

T1 . 
sun 

— 

IT1.  -  T1  | 

1  sun  an 

10 

0.5 

0.3979 

0.4 

0.0021 

10 

0.8 

0.9867 

1.0 

0.0133 

10 

0.9 

1.8734 

2.0 

0.1266 

5 

0.5 

0.2027 

0.2 

0.0027 

5 

0.8 

0.4735 

0.5 

0.0265 

5 

0.9 

0.9851 

1.0 

0.0149 

TABLE  3-3 


T-4475 


PRE  -  T1 

POST  - 

r1 

INITIAL 

p  (p=50) 

APPROX . 

■Hg 

APPROX. 

QUEUE  LENGTH 

EC 

SB 

MEAN=S.D. 

Hi 

1 

HEAN=S.D. 

10 

0.5 

0.0198 

0.0194 

0.02 

0.0405 

0.0405 

0.04 

10 

0.8 

0 . 0204 

0.0206 

0.02 

0.0254 

0.0256 

0.025 

10 

0.9 

0.0201 

0.0200 

0.02 

0.0230 

0.0228 

0.0222 

5 

0.5 

0.0201 

0.0196 

0.02 

0.0405 

0.0406 

0.04 

5 

0.8 

0.0200 

0.0206 

0.02 

0.0255 

0.0254 

0.025 

5 

0.9 

0.0202 

0.0202 

0.02 

0.0229 

0.0228 

0.0222 

*d  =  wan  interdeparture  tiw;  a  atandard  deviation  of  interdeparture  tlw. 


increasing  confidence  to  the  accuracy  and  robustness  of  the 
approximation . 
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3.4 


COMPARISON  WITH  MOORE'S  TECHNIQUE 


In  this  section,  we  provide  some  comparative  results 
between  the  methods  presented  in  this  paper  and  those  proposed 
by  S.C.  Moore  (MOOR  1972]  [MOOR  1975].  As  was  mentioned  in 

Y 

the  introduction,  Moore's  technique  centers  about  the  M  <t)/Ey/l 
queue  and  employs  an  imbedded  Markov  chain  approach  to  analyze 
its  behavior.  In  applying  this  approach,  the  system  status  is 
observed  only  at  departure  instants.  If  a  simple  imbedded 
Markov  chain  approach  (which  assumed  finite  buffer  size)  were 
used,  computational  results  could  be  obtained  very  quickly; 
however,  all  connection  with  the  real  or  continuous  time  axis 
would  be  lost. 

Realizing  this,  Moore  developed  an  exact  technique 
for  maintaining  the  time  axis  connection  in  the  case  of  a  sta¬ 
tionary  arrival  process,  and  an  approximate  technique  for  doing 
the  same  when  the  arrival  process  is  nonstationary.  The  latter 
approximation  resides  in  determining  the  times  of  the  nth  ser¬ 
vices  initiation  and  completion.  If  one  begins  with  Tq  =  0 
and  a  known  queue  length  distribution,  it  is  necessary  to  apply 
the  equations 


+  X 

n 


T  =  t  +  s 

n  n  n 


alternately  in  order  to  maintain  exact  correspondence  between 
the  imbedded  Markov  chain  behavior  (associated  with  the  queue) 
and  the  continuous  time  axis.  In  the  equations,  Tn  is  the 
time  of  the  nth  customers  departure,  Xn  is  the  server  idle 
time  in  the  interval  (T  ^ ,  Tn ) ,  and  Sn  is  the  length  of  the 
nth  customer's  service  time. 
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In  order  to  avoid  the  complexity  of  obtaining  exact 
distributions  for  all  terms  involved,  Moore  develops  approxima¬ 
tions  based  on  the  expected  values  of  and  T  .  These  are 

n  n 

then  used  to  derive  both  the  time  from  departure  n  to  the  next 

arrival  (which  depends  on  T  . )  and  the  number  of  arrivals 

n*  J-  i 

during  service  n  (which  depends  on  Tn>.  The  approximation  is 
due  to  the  fact  that  once  the  arrival  parameters  shift,  the 
expected  values  E[T^]  and  E [ Tfi ]  are  only  approximate. 

Figure  3-4  shows  how  Moore’s  technique  compares  with 
that  developed  in  this  paper.  As  in  Figure  3-3,  the  analysis 
proceeds  over  three  intervals,  the  first  two  of  which  of  2 
second  duration  each  and  the  third  which  is  infinite.  The  p 
associated  with  the  intervals  is  .7,  .8  and  .5,  respectively. 

As  is  seen,  the  two  methods  exhibit  quite  close  agreement. 


4.  SUMMARY  OF  RESULTS  AND  FUTURE  WORK 


In  the  previous  sections,  we  have  examined  in  some 
detail  Stem's  M/M/l/K  approximation  technique  for  obtaining 
the  transient  behavior  of  an  M/M/1  queue  and  proposed  methods 
by  which  it  can  be  applied  in  a  practical,  computationally 
efficient  manner. 

We  have  also  shown  how  the  latter  technique  can  be 
extended  to  accommodate  the  analysis  of  the  M(t)/M/1  queue 
where  we  assume  that  the  M(t)  arrival  process  can  be  approxi¬ 
mated  by  a  discretely  time-varying  Poisson  process.  The  anal¬ 
ysis  of  the  M(t)/M/1  queue  then  becomes  the  analysis  of  a  se¬ 
quence  of  stationary  M/M/l/K  queues  over  finite  intervals  whose 
boundaries  are  defined  by  the  instants  when  the  arrival  param¬ 
eter  changes. 

Together  with  the  above  results,  we  have  proposed  an 
approximate  expression  for  the  departure  process  of  the  M/M/1 
queue  which  implies  that  for  an  M(t)/M/1  queue  whose  arrival 
process  is  discretely  time-varying  M(t),  so  too  the  departure 
process  can  be  approximated  as  discretely  time-varying  (albeit 
of  a  different  form).  The  simulation  results  of  Section  3.3 
have  shown  this  approximation  to  be  accurate. 

Finally,  we  have  compared  our  technique  with  that 
proposed  by  S.C.  Moore  using  an  imbedded  Markov  chain  approach. 
Although  Moore's  technique  permits  the  modeling  of  a  wider 
range  of  queues  than  our  technique  at  this  time,  in  the  partic¬ 
ular  case  of  the  M(t)/M/1  queue  it  is  felt  that  our  technique 
is  more  accurate.  This  is  true  especially  during  transitions 
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to  and  from  intervals  for  which  psl  or  p>l  when  the  arrival 
process  may  change  and  a  number  of  arrivals  occur  between  two 
consecutive  departures.  The  fact  that  our  technique  is  tightly 
coupled  to  the  instants  when  the  parameters  associated  with 
the  underlying  Poisson  processes  change  allows  it  to  "react” 
more  rapidly  to  varying  queue  dynamics. 


The  results  of  this  paper  will  serve  as  the  basis  for 
further  research  into  a  number  of  areas  including: 


•  Extension  and  refinement  of  the  finite 
buffer  approximation  technique  to  more 
general  queues  including  the  M(t)/Ey/1 

Mx(t)/Ey/1  and  E*(t)/Ey/1  queues. 

•  Extension  of  the  methods  and  results  de¬ 
scribed  herein  to  methods  for  analyzing 
feed  forward  and  general  queueing  networks. 

•  Examination  of  the  effects  that  different 
service  disciplines  and  multiple  classes 
have  on  the  behavior  of  such  queues. 

•  Identification  of  the  critical  error 
sources  associated  with  the  approximation 
techniques  and  better  quantification  of 
their  effects. 

•  Application  of  the  methods  to  several 
practical  problems  including  internet¬ 
working  problems,  network  stability  anal¬ 
ysis,  design  and  modeling  of  dynamic  and 
other  routing  strategies,  and  optimal  re¬ 
source  sharing  in  a  transient  environment. 
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